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Abstract 


A  time-dependent  extension  of  the  first-order  hyperbolic  system  method  [J. 
Comput.  Phys.,  227  (2007)  315-352]  for  advection-diffusion  problems  is  in¬ 
troduced.  Diffusive/viscous  terms  are  written  and  discretized  as  a  hyperbolic 
system,  which  recovers  the  original  equation  in  the  steady  state.  The  resulting 
scheme  offers  advantages  over  traditional  schemes:  a  dramatic  simplihcation  in 
the  discretization,  high-order  accuracy  in  the  solution  gradients,  and  orders-of- 
magnitude  convergence  acceleration.  The  hyperbolic  advection-diffusion  sys¬ 
tem  is  discretized  by  the  second-order  upwind  residual-distribution  scheme  in 
a  unihed  manner,  and  the  system  of  implicit-residual-equations  is  solved  by 
Newton’s  method  over  every  physical  time  step.  The  numerical  results  are 
presented  for  linear  and  nonlinear  advection-diffusion  problems,  demonstrat¬ 
ing  solutions  and  gradients  produced  to  the  same  order  of  accuracy,  with  rapid 
convergence  over  each  physical  time  step,  typically  less  than  hve  Newton  iter¬ 
ations. 
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1  Introduction 


In  this  paper,  we  introduce  an  unsteady  hyperbolic  advection-diffusion  scheme. 
The  hrst-order  hyperbolic  system  method  was  initially  proposed  in  Ref.  [1] 
as  a  radical  approach  to  steady  diffusion  problems:  discretize  an  equivalent 
hrst-order  hyperbolic  system  to  solve  the  diffusion  equation.  It  has  a  number 
of  attractive  features  such  as  the  discretization  of  high-order  derivatives  by 
methods  for  hyperbolic  systems,  equal  order  of  accuracy  for  the  solution  and 
the  gradients  (viscous/heat  huxes),  and  orders  of  magnitude  acceleration  in 
steady  convergence.  These  advantages  have  been  demonstrated  for  the  diffu¬ 
sion  equation  [1]  and  the  advection-diffusion  equation  [2]  by  the  second-order 
residual-distribution  method,  and  for  the  compressible  Navier-Stokes  equa¬ 
tions  [3]  by  the  second-order  hnite-volume  method,  and  for  the  advection- 
diffusion  equation  by  hrst,  second  and  third  order  hnite-volume  methods  [4,5]. 
These  lines  of  development,  however,  have  been  restricted  exclusively  to  steady 
state  problems.  Towards  enabling  accurate  practical  time-dependent  computa¬ 
tions  by  the  hrst-order  hyperbolic  system  method,  this  paper  presents  the  hrst 
study  on  the  unsteady  extension  of  the  hrst-order  hyperbolic  system  method. 
We  demonstrate  that  time-accurate  schemes  constructed  via  the  hyperbolic 
method  produce  high-order  accurate  solutions  and  derivatives  at  every  time 
step  and  allow  the  construction  of  a  highly  efficient  inner  solver  for  implicit 
time  stepping. 

Time-dependent  computations  are  possible  by  implicit  time-integration 
methods  such  as  the  second-order  backward-diherence  scheme,  which  is  widely 
used  for  practical  computations.  To  perform  the  implicit  time  integration,  it 
is  required  to  solve  a  system  of  global  time-dependent  residual  equations  over 
each  time  step.  We  employ  a  steady  solver  developed  by  the  hrst-order  hy¬ 
perbolic  system  method  to  solve  the  system  efficiently  and  accurately.  The 
spatial  discretization  is  performed  by  the  residual-distribution  method  as  in 
Refs.  [1,2].  Taking  advantage  of  the  compactness  of  the  residual-distribution 
schemes,  we  develop  a  fully-implicit  steady  solver  with  exact  linearization,  i.e., 
Newton’s  method,  for  a  second-order  upwind  scheme.  The  strategy  taken  here 
may  be  thought  of  as  a  dual-time  formulation  [6].  However,  we  introduce 
a  pseudo-time  just  for  convenience  of  discretization:  the  hrst-order  system 
is  hyperbolic  in  the  pseudo-time  and  thus  will  be  discretized  by  the  upwind 
residual-distribution  scheme.  Once  the  spatial  discretization  is  complete,  the 
pseudo-time  step  is  taken  to  be  inhnity  to  dehne  the  system  of  time-dependent 
residual  equations,  which  is  then  solved  by  Newton’s  method.  In  this  paper, 
we  consider  one-dimensional  linear  and  nonlinear  advection-dihusion  equa¬ 
tions.  These  examples  are  sufficient  to  illustrate  the  general  methodology 
for  enabling  time-dependent  computations  by  the  hyperbolic  method;  this  is 
the  main  goal  of  this  paper.  The  introduction  of  this  unsteady  hyperbolic 
advection-dihusion  scheme  along  with  its  great  advantages  is  served  as  a  foun- 
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dation  for  the  extension  of  the  scheme  to  higher-dimensions  and  more  complex 
systems. 

It  is  emphasized  also  that  there  are  real  practical  applications  that  one¬ 
dimensional  analyses  are  sufficient  and  routinely  being  used  in  industries  (e.g., 
NASA).  For  instance,  in-depth  material  thermal  response  calculations  of  ther¬ 
mal  protection  systems  of  atmospheric  entry  vehicles,  where  advection-diffusion 
equation  is  solved  through  the  material,  is  a  prime  example  of  such  applica¬ 
tions  [7-9].  Another  example  is  the  experimental  aeroheating  data  reduction 
that  are  obtained  with  global  phosphor  thermography  techniques  [10,11].  Ef- 
hcient  and  accurate  one-dimensional  schemes  as  presented  in  this  paper  will 
signihcantly  improve  these  analyses  and  could  potentially  make  such  calcula¬ 
tions  part  of  the  routine  aerothermodynamic  environment  predictions.  The 
present  paper  serves,  therefore,  also  as  an  important  foundation  for  the  exten¬ 
sion  to  more  practical  one-dimensional  problems. 

There  exist  other  methods  for  enabling  time-accurate  computations  by 
the  residual-distribution  method:  the  Crank-Nicolson  method  [12],  the  space- 
time  method  [13-15],  the  backward-difference  method  [16,17],  and  the  explicit 
Runge-Kutta-type  method  [18].  These  methods  are  equally  applicable  to  the 
construction  of  time-accurate  schemes  in  the  hyperbolic  method.  In  this  pa¬ 
per,  we  employ  the  backward-difference  method  for  simplicity  but  without  los¬ 
ing  practical  applicability.  The  main  difference  between  our  time-integration 
method  and  those  in  Refs.  [16, 17]  lies  in  the  method  for  solving  the  sys¬ 
tem  of  residual  equations  arising  from  the  implicit  time  integration  scheme: 
we  employ  Newton’s  method  whereas  Ref.  [16]  employs  the  dual-time  step¬ 
ping  method  with  a  pseudo-time  marching  iteration  and  Ref.  [17]  employs  the 
multigrid  method  with  a  point-implicit  iteration. 

The  paper  is  organized  as  follows.  In  the  next  section,  the  hyperbolic 
advection-diffusion  system  is  described.  In  Section  3,  a  compact  second-order 
residual-distribution  scheme,  a  steady  solver,  and  the  boundary  conditions  are 
discussed.  In  Section  4,  the  construction  of  time-accurate  schemes  is  given. 
In  Section  5,  the  time-accurate  scheme  is  extended  to  a  nonlinear  advection- 
diffusion  equation.  In  Section  6,  numerical  results  are  presented.  Finally, 
Section  7  concludes  the  study  with  remarks. 


2  Hyperbolic  Advection-DifFusion  System 

Consider  the  one-dimensional  (1-D)  advection-diffusion  equation, 

dtu  +  a  dxU  =  V  dxxU,  (1) 

where  a  and  u  are  both  taken  to  be  positive  constant.  We  will  follow  the 
procedure  described  in  Ref.  [2]  and  rewrite  the  above  equation  with  a  hrst- 
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order  advection-diffusion  system: 


(2) 

(3) 


dtu  =  —a  dxU  +  u  dxP, 

dtP  =  {d^u-p)/Tr, 

where  the  relaxation  time,  Tr  >  0,  is  arbitrary.  Towards  the  steady  state, 
the  variable  p  approaches  the  solution  gradient  and  hence  the  above  equation 
becomes  identical  to  the  original  advection-diffusion  equation  in  the  steady 
state.  In  the  vector  form,  the  first-order  advection-diffusion  system  can  be 
written  as 

dtV  +  Ad,V  =  S,  (4) 

where 


a  —v 
-l/Tr  0  J  ’ 


S 


0 

-P/Tr 


(5) 


The  above  system  is  hyperbolic  in  time  because  it  has  the  following  real  eigen¬ 
values  for  any  positive  Tp. 


with  linearly  independent  eigenvectors  [1]. 

Note  that  the  above  hyperbolic  formulation  is  equivalent  to  the  original 
Eq.  (1)  only  in  the  steady  state.  The  time  derivatives  in  the  hyperbolic  sys¬ 
tem,  at  least  should  therefore  be  taken  as  pseudo-time  derivatives.  They 
serve  mainly  as  a  guide  for  discretization  by  making  the  whole  system  hy¬ 
perbolic  in  time  for  which  various  discretization  techniques  are  available,  e.g., 
upwinding.  The  benehts  are,  however,  much  more  than  just  the  convenience  in 
discretization.  First,  the  hyperbolic  discretization  typically  results  in  a  strong 
coupling  between  the  two  variables  and  achieves  the  equal  order  of  accuracy 
for  u  and  p  on  arbitrary  grids.  Second,  the  resulting  explicit  numerical  scheme 
is  stable  with  an  0{h)  time  step  through  the  diffusion  limit  or  0(1/ h)  condi¬ 
tion  number  in  the  residual  Jacobian  for  implicit  solvers.  It  has  been  shown 
to  yield  0(1/ h)  acceleration  in  convergence  over  traditional  methods  for  the 
diffusion  equation  [1,4],  for  the  advection-diffusion  equation  [2,5],  and  for  the 
compressible  Navier-Stokes  equations  [3].  Unsteady  computations  are  possible 
by  implicit  time  stepping,  which  is  the  main  subject  of  this  paper. 

In  the  next  section,  we  hrst  fully  describe  the  advection-diffusion  equation 
discretization  process,  the  implicit  steady  state  formulation,  and  the  boundary 
condition  implementation  in  the  steady  state  limit.  We  then  extend  the  scheme 
to  time-accurate  and  nonlinear  advection-diffusion  equation  in  the  sections 
that  follow. 
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Figure  1:  Schematic  of  grid  spacing  for  a  1-D  grid. 


3  Discretization,  Steady  Solver,  and  Bound¬ 
ary  Condition 

3.1  Discretization 


To  discretize  the  hrst-order  hyperbolic  advection-diffusion  system,  we  employ 
Residual-Distribution  (RD)  method.  The  method  consists  of  two  steps;  1) 
residual  evaluation  over  the  cells  (or  elements),  and  2)  distribution  of  the 
residuals  to  the  nodes. 

Consider  a  one-dimensional  domain  discretized  with  N  nodes  that  are  dis¬ 
tributed  arbitrarily  over  the  domain  of  interest  with  the  solution,  u,  and  the 
solution  gradient,  p,  data  stored  at  the  nodes  denoted  by  Xi,  i  =  1,2,  3, ...,  N 
(Fig.l).  We  define  the  cell-residual  by  integrating  the  spatial  part  of  the 
Eq.  (4)  over  the  cell,  C,  defined  by  the  nodes  i  and  i  +  1\ 


fXi  +  l 


-a{ui+i  -  Ui)  z/(pi+i  -  Pi) 


= 


-A^xU  -I-  S)dx  = 


T 

±  r 


Pi+l  +  Pi  ,  N 


(7) 

Note  that  the  source  term  integration  has  been  evaluated  by  the  trapezoidal 
rule  to  ensure  second-order  accuracy.  It  should  be  noted  also  that  the  flux  bal¬ 
ance  terms  e.g.,  u^,  has  been  integrated  exactly,  which  is  not  possible  in  higher 
dimensions  and  a  high-order  quadrature  would  be  required  to  achieve  the  equal 
order  of  accuracy  for  the  solution  and  the  gradients  in  the  residual-distribution 
method  as  described  in  Ref.  [2].  We  again  remark  that  the  relaxation  time, 
Tr,  is  arbitrary  because  in  the  steady  state  limit  p  =  for  any  T^.  We  com¬ 
plete  the  evaluation  of  the  cell  residual  with  the  definition  of  T^.  Using  the 
dimensional  analysis  [2],  this  is  defined  as 


T 

J.  rjr. 


max(Ai,  A2) 


a 

2 


(8) 


where  Lj.  is  a  length  scale  that  may  be  optimized  using  Fourier  analysis  to 
enhance  the  convergence  [19].  The  effect  of  the  optimized  length  scale  remains 
to  be  investigated.  Here  we  use  the  value  recommended  in  Ref.  [1,4,5]: 


Pj  rp  - 


1 


(9) 
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We  now  solve  for  from  Eq.  (8): 

T  = 

-L  r 


(10) 


a  +  z//Lr 

Substituting  back  the  relaxation  time  into  Eq.  (6),  we  express  the  eigenvalues 


as 


Xi  =  —h'/Lr,  \2  =  a  +  v/Lr.  (11) 

At  this  point  the  residual  evaluation  is  completed.  We  now  dehne  the 
distribution  of  the  residuals.  We  achieve  this  by  first  diagonalizing  the  matrix 
A  with  the  following  right  and  left  eigenvectors: 


R  = 


— AiT^  —X^Tj. 

^  (xLy  ~\-  y 

1  -Lr 

1  1 

—  1  {uLr)  /  {aLr  +  v) 

,(12) 


We  then  expand  the  matrix  A  to  two  separate  parts  that  are  distinguished  by 
its  corresponding  wave  speeds  (eigenvalues): 


where 


A  =  RAL  =  XiR-  +  A25+, 


A 


Ai  0 
0  A2  ’ 


(13) 

(14) 


R- 


1 

QjL'p  ~\-  2z/ 


z/  z/L^ 

{ojL^  ~\-  j QjL'p  -\-  y 


(15) 


(iLj>  ‘2iy 


aLr  +  z/  —uLr 
—  {aLr  +  l')/Lr  V 


(16) 


The  matrices  and  B~  project  the  solution,  respectively,  onto  the  left- 
and  right-running  waves.  Hence,  this  is  upwinding.  We  also  note  that  the 
distribution  matrices  have  the  following  property 


=  /, 


(17) 


which  is  required  for  conservation.  We  can  now  distribute  the  cell  residual 
to  the  nodes  using  the  projection  matrices  B^  and  B~  as  described  in 
Fig.  2.  The  distribution  is  done  this  way  because  the  left  running  wave,  A2, 
for  example,  is  contributing  to  the  B^  and  thus  the  cell  residual  on  the  left  cell 
is  weighted  with  the  B^ .  After  the  distribution  step,  we  obtain  the  following 
semi-discrete  scheme  in  pseudo-time: 

^  =  1(5+$^  +  H-$^),  (18) 

dr  hi 

where  and  denote  the  cell-residuals  over  the  left  and  right  cells  of  the 
node  i,  respectively,  and  hi  is  the  dual  volume  (see  Fig.  1)  defined  by 


hi 


hi  +  hji 


hj^  Xi  /i/j  Xi- 


(19) 
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Figure  2:  Residual  distribution  to  the  nodes. 


3.2  Steady  Solver 

Equation  (18)  can  be  solved  with  either  an  explicit  or  an  implicit  formula¬ 
tion.  We  choose  to  focus  on  the  implicit  formulation  for  a  better  transition 
to  the  next  sections,  where  we  extend  this  scheme  to  general  nonlinear  time- 
dependent  advection-diffusion  problems.  Therefore,  we  drop  the  pseudo-time 
derivative  from  (18)  and  solve  the  resulting  system  of  steady  residual  equa¬ 
tions. 


0  =  Res, 


(20) 


by  a  fully  implicit  solver.  We  remark  that  in  some  complex  problems  the 
pseudo  time  derivative  should  be  kept  for  stability  purposes.  In  our  one¬ 
dimensional  problems,  however,  we  could  choose  an  inhnite  pseudo  time  step 
without  encountering  any  stability  issue. 

The  implicit  formulation  is  dehned  by 

U^+i  =  +  AU^  (21) 

where  U  =  {ui,pi,U2,P2,  ■  ■  ■  ,un,Pn)  and  k  is  the  iteration  counter.  The 
correction  AU^  =  is  determined  as  the  solution  to  the  linear 

system: 


=  -Res‘,  (22) 

where  Res^  is  the  steady  residual  vector  evaluated  by  U^.  The  Jacobian 
matrix  is  sparse,  having  three  2x2  blocks  in  each  row  for  all  interior  nodes 
and  two  blocks  for  boundary  nodes.  The  Tth  interior  pair  of  row  of  the  linear 
system  is  given  by 


+  J.Af/f  +  (23) 

fti 


where  =  (Awf,  Apf),  A?7f+i  =  (Awf+i,  Apf+i), 


Ji—l 

J^ 

Ji+1 


d<l>^ 

dU~i^ 


1 

hi 


(9$^ 

m 


+  B- 


~m)  ’ 


1  9$^ 

hi  dU ij^\ 


and  the  derivatives  of  the  cell-residuals  are  given  by 


(24) 

(25) 

(26) 

(27) 


r  —a  V 

Wl  ~  ^ 

L  Tr  2Tr 


(28) 


The  linear  system  may  be  solved  efficiently  by  Thomas’  algorithm,  which 
is  an  0{N)  method.  In  this  paper,  however,  we  employ  the  Gauss-Seidel  (GS) 
relaxation,  which  is  also  an  0{N)  method  for  the  discretization  arising  from 
hyperbolic  advection-diffusion  system  and  extends  straightforwardly  to  more 
complex  systems,  irregular  grids,  and  higher  dimensions,  whereas  Thomas’ 
algorithm  is  not  applicable  in  higher  dimensions.  The  GS  relaxation  scheme 
may  be  applied  equation  by  equation  (decoupled  relaxation)  or  node  by  node 
(collective  relaxation)  updating  the  solution  set  at  a  node  simultaneouously 
[20].  In  this  work,  the  collective  GS  relaxation  is  employed  for  robustness;  the 
decoupled  GS  relaxation  was  found  to  give  better  convergence  in  most  cases 
but  also  found  to  be  unstable  in  some  cases  whereas  the  collective  relaxation 
encountered  no  such  problems.  In  actual  computations,  we  do  not  fully  solve 
but  relax  the  linear  system  to  reduce  the  linear  residual  by  two  or  three  orders 
of  magnitude. 

Note  that  the  Jacobian  is  exact  and  thus  the  implicit  algorithm  is  New¬ 
ton’s  method.  It  is  one  of  the  advantages  of  the  residual-distribution  method 
that  second-order  accuracy  is  achieved  within  a  compact  stencil,  thus  allowing 
the  exact  linearization  with  a  sparse  Jacobian  matrix.  The  same  is  true  for 
extension  to  higher  dimensions.  For  linear  problems,  the  method  is  essentially 
a  direct  method.  That  is,  if  we  solve  the  linear  system  exactly,  then  we  get 
the  solution  at  the  hrst  Newton  iteration. 
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Figure  3:  Schematic  of  the  Dirichlet  and  the  Neumann  boundary  conditions 
implementation  in  the  hyperbolic  advection-diffusion  scheme. 


3.3  Boundary  Condition 

In  this  section  we  describe  the  boundary  condition  implementation  for  our  hrst- 
order  hyperbolic  system  scheme.  Because  in  our  hyperbolic  advection-diffusion 
system  the  second  equation  in  the  hyperbolic  system  solves  for  the  solution 
gradient,  p,  the  Neumann  boundary  condition  is  implemented  exactly  in  the 
same  way  as  the  Dirichlet  boundary  condition.  Therefore,  we  only  show  the 
formulation  for  the  u  variable  and  the  same  can  be  repeated  for  the  p  variable. 
Note  that  the  discrete  problem  has  a  unique  solution  for  linear  problems  with 
two  values  given  at  the  left  and  right  boundaries  (this  is  consistent  with  the 
differential  equation  allowing  two  boundary  conditions):  2(iV— 1)  cell-residuals 
for  2N  —  2  unknowns,  which  is  the  case  in  both  of  the  following  boundary 
conditions. 


Dirichlet  or  Neumann  Boundary  Condition 


For  Dirichlet  or  Neumann  boundary  condition  type  schematically  shown  in 
Fig.  3,  the  discretized  equation  for  the  boundary  condition  imposed  on  u  takes 
the  following  form  at  the  hrst  node: 

ADf  +  J;AU^  =  (29) 

hi 

where  for  the  imposed  boundary  condition  on  the  u  variable,  the  Jacobian 
matrices  are 


J2 


1  0 

d<h^(2)  d$^(2)  ’ 

dui  dpi 

0  0 

d$^(2)  d<h^(2)  ’ 

du2  dp2 


(30) 


(31) 


where  <F'^(2)  denotes  the  second  component  of  corresponding  to  the  second 
equation  (the  p  equation)  of  the  hrst-order  hyperbolic  system.  The  same 
process  is  repeated  for  imposed  boundary  condition  on  the  last  node  (for  either 
u  or  p.) 
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Figure  4:  Schematic  of  the  periodic  boundary  condition  implementation  in  the 
hyperbolic  advection-diffusion  scheme. 


Periodic  Boundary  Condition 

To  implement  a  periodic  boundary  condition  we  set  Utv  =  Ui  in  Eq.  (23)  and 
obtain  the  following  discretized  equation  and  applied  it  on  both  two  boundaries 
(Fig.  4): 


Jn-iAU^n-i  +  Ji^U^  +  J2^U^  =  (32) 

til 

where  L  denotes  the  cell  dehned  by  the  nodes  iV  —  1  and  N,  and  R  denotes 
the  cell  dehned  by  the  nodes  1  and  2. 


4  Extension  to  Time-Dependent  Problems 

We  extend  the  steady  state  formulation  by  modifying  the  hrst-order  hyperbolic 
system  given  by  Eq.  (2)  and  rewriting  it  as  a  dual-time  step  system: 

d-rU  =  —aUx  -l-  i^dxP  —  -^u  +  s,  (33) 

drP  =  idxU-p)/Tr,  (34) 


where  t  is  the  physical  time,  At  is  the  physical  time  step,  r  is  the  pseudo-time, 
and  s  and  a  both  depend  on  the  physical  time-step  discretization  scheme. 
Adopting  hrst-  or  second-order  Backward-Differencing-Formulation  (BDFl  or 
BDF2,  respectively)  for  the  physical  time  integration,  the  a  and  s{x)  terms 
are: 

a  =  1,  s  =  —  :  BDFl, 


_  2AC-^  +  AC  _  AP-^  +  AC  , 

~  +  Ai"  ’  ^  Ai’*(Ai"-i  +  Ai")^  '  ^ 

(35) 

where  a  nonuniform  time  step  formulation  is  used  in  the  BDF2  formulation. 
The  BDFl  is  only  used  at  the  very  beginning  of  the  simulation  with  an  initial 
time  interval  of  At^{<<  AD).  After  the  hrst  time  interval,  BDF2  with  a 
uniform  time  step.  At,  is  used.  This  procedure  will  ensure  a  second-order 
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accurate  result  through  all  times.  The  use  of  BDFl  is  necessary  because 
explicit  time  stepping  is  not  possible  with  the  hyperbolic  system  method. 

We  solve  the  above  system  over  each  physical  time  interval  to  obtain  the 
steady  state  solution  in  the  pseudo-time  space.  In  the  steady  state  in  r,  we 
have  p  =  dxU  and  therefore  recover  the  consistent  implicit  time  discretization 
scheme  for  the  original  time-dependent  advection-diffusion  equation. 

In  the  dual-time  step  system,  the  physical  time  derivative  term  acts  as  a 
source  term  to  the  hrst  equation.  This  system  thus  shares  the  same  eigenvalues 
and  eigenvectors  as  the  steady  hyperbolic  advection-diffusion  system.  We  will 
see  this  more  clearly  by  writing  the  time-dependent  hyperbolic  system  in  the 
vector  form; 


(9U  9U 

dr  dx 


DU  +  S, 


(36) 


where 


A 


r  a  —v 

,  D  = 

0 

c  _ 

s{x) 

1 - 

1 

o 

0 

-1/T,  _ 

)  ^  — 

0 

(37) 


Because  of  the  similarity  in  wave  dehnitions  in  both  time-dependent  and  steady 
hyperbolic  systems,  the  distribution  matrices  B~  and  are  the  same  for  both 
of  these  systems.  Once  the  discretization  is  obtained,  we  drop  the  pseudo-time 
derivatives  and  solve  the  resulting  system  of  discrete  equations  by  Newton’s 
method  as  described  in  Section  3.2.  To  develop  a  time- dependent  scheme,  we 
therefore  only  need  to  add  the  source  term  effect  in  the  cell-residual. 

The  cell-residual  of  the  time-dependent  hyperbolic  advection-diffusion 
system  is 


= 


-AU^  +  DU  +  S)dx 


a 


-a{ui+i  -  Ui)  +  iy{pi+i  -  Pi)  -  {xi+i  -  Xi)  —  {ui+i  +  Ui)/2 


T 

±  r 


Pi+i  +  Pi  I 

'^i+l  2  (^i+1 


k+l 


+ 


{Xi+i  -  Xi){si+i  +  Si)l2 


0 


(38) 


where  k  and  n  are  the  Newton  iteration  counter  and  the  physical  time  index, 
respectively.  Note  that  the  second  term  of  the  Eq.  (38),  which  is  computed  at 
the  two  previous  physical  time  steps,  is  constant  during  the  Newton  iteration, 
and  thus  will  not  contribute  to  the  Jacobian.  The  residual  Jacobians  needed 
in  the  Newton  solver  are  exactly  in  the  same  form  as  Eq.  (26),  but  the  deriva¬ 
tives  of  the  cell-residuals  now  include  the  contribution  from  the  physical  time 
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derivative: 


a 


—  V 


dU, 


d^L 

dUi 


a  —  Hr 


2At 


1 

±  r 

—a  —  hi 


a 


1 

±  r 


hR  ’ 
2Tr  . 

V 

. 

'Wr  . 


(39) 


(40) 


The  hyperbolic  advection-diffusion  scheme  is  now  time-accurate.  At  each 
physical  time  level  n,  we  solve  the  pseudo  steady  problem  by  Newton’s  method 
with  the  current  solution  at  n  as  the  initial  solution. 


5  Extension  to  Nonlinear  Advection-Diffusion 
System 

Consider  the  generalized  time-dependent  nonlinear  advection-diffusion  equa¬ 
tion: 

dtu  + d^f  =  (41) 

where  /  is  a  nonlinear  function  of  u  and  v  =  z/(m)  is  a  function  of  u.  We  write 
the  above  equation  using  the  dual-time  stepping  formulation  in  the  following 


first-order  hyperbolic  system  as 

drU  =  -d^f  A  dxP  -  -^u  +  s, 

(42) 

Tr 

,  .drP  =  d^u  p/v{u). 
v{u) 

(43) 

This  is  a  preconditioned  conservative  form  introduced  in  Ref.  [3]  to  extend  the 
hyperbolic  method  to  nonlinear  equations.  Note  that  the  variable  p  will  be 

the  diffusive  flux  in  the  pseudo  steady  state,  not  the  gradient, 
form,  the  preconditioned  conservative  system  is  written  as 

In  the  vector 

P  iU.  =  -F,  +  S, 

(44) 

where 

P  ^ 


1  0 
0  Tr/v{u) 


F 


f -p\ 

-u 


S 


-^u  +  s{x)  \ 

-v/y{u)  )  ’ 


(45) 


and  P  is  the  preconditioning  matrix.  The  main  role  of  the  preconditioning 
matrix  here  is  not  to  optimize  the  condition  number  of  the  differential  system, 
but  rather  to  simplify  the  discretization.  In  particular,  in  this  form,  there  is 
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no  need  to  differentiate  the  diffusion  coefficient  with  respect  to  the  solution  to 
derive  the  flux  Jacobian.  Multiplying  both  sides  by  P  from  the  left,  the  flux 
term  becomes 

OF 

.  (46) 


PFx  =  Pt^Ux  =  AUx 


where 


A  = 


(9U 

a 

-yiT, 


-1 
^  0 


(47) 


Here,  we  have  introduced  a  =  Of  /du,  which  is  not  a  constant  but  a  func¬ 
tion  of  u.  Note  that  the  preconditioned  Jacobian  is  formally  equivalent  to 
the  Jacobian  of  the  linear  time  dependent  hyperbolic  advection-diffusion  sys¬ 
tem.  Hence,  the  eigenvalues  and  eigenvectors  of  this  system  are  identical  to 
those  of  the  linear  equations.  It  greatly  simplihes  the  construction  of  the  dis¬ 
tribution  matrices  as  described  below.  This  is  one  of  the  advantages  of  the 
preconditioned  conservative  formulation. 

To  dehne  the  cell-residual,  we  hrst  integrate  the  right  hand  side  of  Eq.  (44) 
to  get  the  unpreconditioned  residual  T'"; 


= 


r^i+i 


(  —  Fa;  -f  S)dx 


'  Xi 


=  -(Fa+i-Fa)  +  ^^±^^(a;,+i-a;a).  (48) 

For  the  advective  flux  that  is  quadratic  in  the  solution,  it  is  well  known  that 
a  conservative  linearization  exists: 


+  (49) 

where  is  the  analytical  Jacobian  evaluated  at  the  arithmetic  average  of  the 
solution,  (7/j+i  -|-  Ui)/2  and  it  satisfies 


OF 


Fi+i  ~  Fj  —  ^u(Ui+i  —  Uj), 


(50) 


exactly.  The  preconditioned  cell-residual,  which  is  to  be  distributed,  is  then 
given  by 


c 


(51) 


where  the  matrix  P  is  evaluated  at  the  arithmetic  average  of  the  solution  so 
that  we  have 


=  -A(U,+i-U0  + 


Si+i  J-  Sj 


{Xi+i  -  Xi). 


(52) 
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The  distribution  matrices  can  then  be  constructed  in  the  same  way  as  in  the 
linear  case  based  on  A; 


B7.  = 


C  ~ 


vLr 


(yOjLrf  T  j  (XLrp  -\-  V 


(53) 


1 

(xL'f‘  ~\~  2z^ 


aLr  +  v  —uLr 
—  {aLr+V)/Lr  V 


(54) 


where 


a  = 


O-i  +  O-i+l 
2 


z/  = 


+  ^i+1 
2 


(55) 


The  subscript  C  indicates  that  these  matrices  are  dehned  over  the  cell  C. 
Note  that  the  distribution  matrices  are  not  globally  constant  in  general  for 
nonlinear  equations.  The  construction  of  the  distribution  matrices  based  on 
the  conservative  linearization  is  important  for  proper  upwinding.  If  different 
averages  are  used,  the  distribution  may  not  be  strictly  upwind.  If  the  target 
equation  does  not  allow  the  conservative  linearization,  the  construction  pro¬ 
posed  in  Ref.  [21]  may  be  employed  to  design  a  proper  upwind  distribution 
matrix.  Note  also  that  we  evaluate  u  in  the  source  term  by  the  average  value 
over  the  cell  to  simplify  the  second  component  of  the  cell-residual  as  follows: 


$ 


c 


fi)  +  {pi+1  -  Vi)  -  {Xi+i 


Xi) 


a 

At 


("Wi+l  T 


k+1 


^  [u{Ui+i  -  Ui)  -  p{Xi+i  -  Xi) 


+ 


-f-  Si)  12, 


0 


n— l,n 

,  (56) 


where 


T 

J.  ^ 


a  +  V/  Lp 


(57) 


Observe  that  we  have  p  =  V{ui+i  —  Ui)/{xi+i  —  Xi)  in  the  pseudo  steady  state. 
The  cell-residual  is  then  distributed  to  the  nodes  in  all  cells  to  yield  the  residual 
equation  at  node  v. 

0  =  (58) 

where  the  pseudo-time  derivatives  have  been  dropped. 


15 


We  note  that  the  Jacobian  of  this  nonlinear  advection-diffusion  system 
is  slightly  different  than  that  of  the  linear  problem.  The  Jacobian  blocks, 
however,  are  function  of  derivative  of  the  distribution  matrices  B~  and 


Ji—l 


J^ 


Jr 


i+l 


1  d  (5+$^) 

hr  ’ 

h,  au,  au,  j ' 
1  a 

hi  OUi^\ 


(59) 

(60) 
(61) 


In  the  special  case  ajv  =  constant,  which  is  the  only  case  considered  in  this  pa¬ 
per,  the  distribution  matrices  reduce  entirely  to  constant  matrices.  Therefore, 
we  can  entirely  ignore  the  derivatives  of  B~  and  B'l  in  deriving  the  Jaco¬ 
bian.  As  a  result,  the  Jacobians  are  given  formally  as  in  Eq.  (26).  This  makes 
the  implementation  of  both  linear  and  nonlinear  systems  extremely  easy  and 
straightforward.  Effects  of  nonlinearity  arise  only  in  the  derivatives  of  the  left 
and  right  cell-residuals: 


Ikl 


o-i  —  hn 


Ui+i  -Ui  du  u{ui+i  -  Ui) 


T 

±  r 


dui 


+ 


a 

phR 


-1 


k 


da 


+ 


1  dv 


duj  Lr  duj 


V 


T 

±  r 


h 


R 


2Tr  J 

(62) 


d^^ 


~Ri 


hr 


a 


k 


U, 


‘2  At 

-Ui-i  dV  v{ui-  Ui-i)  -  phi  (  da 


T 

±  r 


du. 


+ 


- 


+ 


1  du 


V  dui  Lr  dui 


V 

+  = 
Tr 


hr 


2Tr  J 

(63) 

where  the  average  values  indicated  by  the  over  bar  are  understood  to  be  taken 
over  the  corresponding  cell.  Note  that  these  derivatives  depend  on  the  solution 
at  the  current  iteration,  and  therefore  the  superscript  k  has  been  introduced. 
The  same  applies  to  the  distribution  matrices  also.  Obviously,  the  nonlinear 
system  approaches  the  linear  system  in  the  event  that  both  the  advection  and 
the  diffusion  coefficients  are  constant.  Also,  the  boundary  conditions  can  be 
implemented  as  described  in  Section  3.3  for  the  linear  case.  Therefore,  our 
hyperbolic  nonlinear  time  dependent  advection-diffusion  system  scheme  can 
be  implemented  for  linear,  nonlinear,  steady  and  unsteady  problems  without 
loss  of  accuracy.  Eor  simplify  the  discussion  we  only  showed  Jacobian  for  non¬ 
linear  cases  with  ajv  =  constant  but  the  scheme  works  for  general  nonlinear 
problems. 
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6  Results 


In  this  section  we  present  the  results  in  three  categories:  1)  Steady  advection- 
diffusion  equation  for  high  Reynolds  (or  Peclet)  number  applications,  2)  Un¬ 
steady  linear  advection-diffusion,  and  3)  Unsteady  nonlinear  advection-diffusion 
problems.  We  then  present  the  order  of  accuracy  results  in  the  last  subsec¬ 
tion.  Note  that  time-accurate  computations  are  started  by  BDFl  in  the  first 
time  step  with  extremely  small  time  step  (e.g.  t  =  10“®),  and  then  by  BDF2 
thereafter  with  much  larger  time  step  (e.g.  t  =  0.001  —  0.5)  for  all  unsteady 
cases.  This  will  ensure  that  the  order  of  accuracy  stays  second  order  through 
all  times.  We  remark  that  explicit  time  stepping  is  not  available  for  time- 
accurate  computations  with  the  hyperbolic  system  method. 


6.1  Steady  Linear  Advection-Diffusion  (High  Re  Appli¬ 
cations) 

Consider  the  advection-diffusion  equation  in  a;  G  (0, 1)  with  m(0)  =  0  and 
m(1)  =  1  boundary  conditions: 

dtu  +  adxU  =  fdxxU  -|-  s(a;),  (64) 


where 


TT 


s(a;)  =  — [acos(7ra;)  -|- 7rz/sin(7ra;)],  Re  =  ajv. 
Re 


(65) 


This  is  a  boundary  layer  problem  with  a  non-trivial  steady  state  solution  in 
the  diffusion  limit  as  a  result  of  the  source  term  addition  [2]  .  This  equation 
develops  a  very  narrow  boundary  layer  near  the  right  boundary  (a;  =  1)  where 
the  advection  becomes  dominant.  The  exact  steady  state  solution  to  this 
problem  is  given  by 


^—Re 


U 


exact  I 


^{x—l)Re 


X  = 


o—Re 


+  —  sin(7ra;). 
Re 


(66) 


We  chose  Re  values  ranging  from  1  to  10®  and  solved  the  equation  on  nonuni¬ 
form  grid  sizes  up  to  30000  nodes.  We  remark  that  the  high-i?e  cases  re¬ 
quired  extremely  hne  grids  to  meet  the  well-known  requirement  on  the  mesh 
Reynolds-number  [2].  If  desired,  the  computations  can  be  performed  on  sub¬ 
stantially  coarser  grids  with  more  aggressive  and  customized  grid  stretching 
as  well  as  with  some  nonlinear  mechanism  to  prevent  numerical  oscillations. 
However,  we  simply  rehned  the  grid  to  meet  the  mesh  Reynolds-number  re¬ 
quirement  because  our  method  is  powerful  enough  to  solve  the  problem  very 
efficiently  (i.e.,  less  than  10  Newton  iterations)  even  on  such  dense  grids.  The 
ability  to  efficiently  solve  the  problem  on  highly  rehned  grids  is  a  great  ad¬ 
vantage  of  these  schemes.  Shown  in  Fig.  5  are  the  comparisons  between  the 
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Figure  5:  Residual-Distribution  scheme  on  hyperbolic  advection-diffusion  sys¬ 
tem  for  boundary  layer  equation. 


computational  and  the  exact  solutions.  In  this  comparison,  the  exact  analytic 
solution  is  represented  with  -|-  (plus)  symbols  while  the  numerical  results  are 
shown  with  solid  lines  (denoted  Hyp-ADE).  The  exact  solutions  are  plotted  at 
nodes,  and  therefore  they  also  show  the  computational  grids.  The  solutions 
were  obtained  with  the  Newton-GS  method  as  described  in  Section  3.2.  For 
this  problems,  within  each  Newton  iteration,  the  GS  relaxation  were  conducted 
until  three  orders  of  magnitude  reduction  is  achieved  in  the  linear  solver.  The 
computations  were  continued  until  the  residuals  of  both  the  solution  and  the 
solution  gradient  were  reduced  by  eight  orders  of  magnitude.  Table  1  gives 
the  convergence  data  obtained  using  the  implicit  formulation  of  the  hyper¬ 
bolic  advection-diffusion  scheme.  The  linear  dependency  of  the  iterations  on 

Table  1:  Boundary  layer  problem  (Gonvergence  criteria:  Residuals  <  10“®.) 


logio 

Number  of  nodes 

GS  relaxations/Newton  iteration 

Newton  iteration 

0 

10 

17 

7 

0 

100 

324 

7 

0 

200 

631 

7 

0 

300 

967 

7 

1 

300 

549 

7 

2 

300 

131 

7 

3 

300 

28 

7 

4 

300 

38 

6 

5 

3000 

60 

6 

6 

30000 

50 

5 

the  grid  size  is  demonstrated  at  Re  =  1,  which  is  a  consequence  of  solving  the 
advection-diffusion  equation  as  a  hyperbolic  system  for  all  Reynolds  numbers 
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through  the  diffusion  limit.  This  is  remarkable  because  the  linear  convergence 
is  retained  for  any  irregular  grid  in  any  dimensions  {N  is  approximately  the 
number  of  nodes  in  each  coordinate  direction  in  two  and  three  dimensions) 
as  demonstrated  in  Refs.  [1-5].  In  conventional  methods,  e.g.,  such  as  the 
alternating  direction  implicit  method  with  Thomas’  algorithm,  the  linear  con¬ 
vergence  can  be  achieved  only  on  Cartesian  grids  (one-dimensional  grid  must 
exist  in  each  coordinate  direction).  The  0{N)  convergence  on  general  prob¬ 
lems  and  arbitrary  grid  systems  makes  the  current  scheme  extremely  attractive 
because  it  leads  to  orders  of  magnitude  faster  convergence  in  comparison  with 
conventional  methods  whose  convergence  is  typically  O(iV^). 

Table  1  shows  for  the  grid  of  300  nodes  that  the  number  of  GS  relaxations 
decreases  signihcantly  as  Re  increases.  In  Ref.  [2] ,  a  formula  for  L^.  was  derived 
that  minimizes  the  condition  number  of  the  hyperbolic  advection-diffusion 
system  and  a  nearly  uniform  convergence  was  obtained  for  an  explicit  pseudo¬ 
time  marching  scheme.  We  tested  the  formula  in  our  method,  but  found  that  it 
did  not  change  the  signihcant  decrease  although  it  did  reduce  the  number  of  GS 
relaxations  approximately  by  10  %.  A  different  approach  may  be  necessary  to 
derive  a  formula  to  achieve  Re-independent  convergence  of  the  GS  relaxation 
in  the  implicit  solver.  Further  study  on  Re-independent  convergence  is  left 
as  future  work.  Nevertheless,  it  should  be  noted  that  the  number  of  Newton 
iterations  is  nearly  Re-independent  and  very  small:  only  5  to  7  iterations  to 
reduce  the  residual  by  eight  orders  of  magnitude.  Gonsidering  the  fact  that  the 
cost  of  one  GS  relaxation  is  signihcantly  cheaper  than  one  Newton  iteration, 
we  hnd  that  the  developed  solver  is  still  a  very  powerful  steady  solver.  It  thus 
serves  as  a  very  efficient  solver  for  the  unsteady  residual  equations  as  shown 
in  the  following  sections. 

Finally,  we  remark  that  the  high-Re  cases  required  extremely  hne  grids  to 
meet  the  well-known  requirement  on  the  mesh  Reynolds-number  as  described 
in  Ref.  [2].  If  desired,  the  computations  can  be  performed  on  substantially 
coarser  grids  with  more  aggressive  grid  stretching.  However,  we  simply  rehned 
the  grid  to  meet  the  mesh  Reynolds-number  requirement  because  our  method 
is  powerful  enough  to  solve  the  problem  very  efficiently  (i.e.,  5  to  7  Newton 
iterations)  even  on  such  dense  grids.  The  ability  to  efficiently  solve  the  problem 
on  highly  rehned  grids  is  a  great  advantage. 


6.2  Unsteady  Linear  Advection-Diffusion 

We  present  unsteady  computations  for  the  two  types  of  boundary  conditions; 
Dirichlet  (or  Neumann)  and  periodic.  We  again  note  that  for  hyperbolic  sys¬ 
tem  scheme  the  implementation  of  Dirichlet  and  Neumann  boundary  condi¬ 
tions  are  similar. 


19 


Consider  the  time-dependent  advection-diffusion  equation  in  a;  G  (0, 1) 


dtu  +  adxU  =  vdxxU,  (67) 

with  the  following  time- dependent  boundary  condition  (Dirichlet): 

M(0,t)  =  0,  (68) 

u{l,t)  =  U  cos{ut),  (69) 


where  U  is  the  amplitude  of  the  oscillation  and  u  is  the  frequency  of  the  os¬ 
cillation  on  the  right  boundary.  This  problem  has  the  following  exact  analytic 
solution: 


u 


exact 


{x,  t)  =  Real 


^Xix 


^\2X 


oX\ 


0X2 


-Ue 


iiot 


a  ±  \/a^  R  Mujv 

’  “  2z^  ’  ^  ^ 


where  i  =  We  solved  this  time-dependent  advection-diffusion  equation 

with  the  following  first-order  hyperbolic  system: 


drU  -|-  adxU  = 

^  3 

VOxP - -I - - , 

^  2At  2At 

(71) 

drP  = 

{dxU-p)/Tr. 

(72) 

The  solution  evolution  was  started  with  the  exact  initial  solution  given  as 

(„Aix  _  „\2X  \ 

— in - 

and  was  advanced  in  physical  time  after  the  time-dependent  residual  was  re¬ 
duced  by  two  orders  of  magnitude  at  each  physical  time  step.  We  tested  the 
scheme  on  a  number  of  nonuniform  grid  systems  ranging  from  iV  =  100...  10000 
nodes.  Shown  in  Fig.  6  are  the  solution  and  the  solution  gradient  obtained 
on  the  coarsest  mesh  (N=100)  for  f/  =  1,  z/  =  a  =  1,  and  uj  =  77r/2.  We 
show  the  results  on  the  coarsest  mesh  just  to  demonstrate  the  accuracy  of  the 
scheme.  Similar  results  were  obtained  on  hner  grid  systems.  The  results  are 
then  compared  with  the  analytic  solution.  Here  the  numerical  results  (denoted 
Hyp-ADE)  are  over  plotted  with  the  exact  solutions,  meaning  excellent  agree¬ 
ment.  We  will  show  the  residual  values  in  the  accuracy  verification  section. 

We  also  obtained  convergence  data  for  different  grid  systems  over  100  time 
steps.  These  are  tabulated  in  Table  2,  which  shows  a  linear  increase  in  number 
of  iterations  with  increase  of  the  grid  system.  This  is,  again,  a  consequence  of 
solving  the  advection-diffusion  equation  as  a  hyperbolic  system,  where  there 
are  only  hrst-derivatives  and  therefore  the  number  of  iterations  is  0{N),  not 
0(7V^).  It  is  also  shown  in  Table  2  that  only  four  Newton  iterations  are 
required  for  all  grids  considered. 


20 


Figure  6:  Time-dependent  advection-diffusion  with  oscillatory  boundary  con¬ 
dition  u{0,t)  =  0,  u{l,t)  =  cos(^t)  on  iV  =  100  irregular  nodes.  Time  step  = 
0.01s.  Left  figures  at  f  =  0.1s,  right  hgures  at  f  =  0.4s. 


Table  2:  Hyperbolic  time-dependent  linear  advection-diffusion  problem  with 
BDF2.  Average  data  over  100  time  steps  are  given.  (Convergence  criteria: 
Residuals  <  10“^) 


Number  of  nodes 

GS  relaxations/Newton  iteration 

Newton  iterations 

100 

127 

4 

300 

370 

4 

500 

607 

4 

1000 

1201 

4 
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For  the  periodic  boundary  condition  demonstration,  consider  the  time- 
dependent  advection-diffusion  equation  in  a:  G  (0,1)  given  as  in  Eq.  (67)  with 
the  following  initial  condition: 

M(a;,  t  =  0)  =  sin(fi;a;),  (74) 

where  k  is  an  arbitrary  constant.  The  exact  solution  to  this  problem  with  a 
periodic  boundary  condition  is: 

t)  =  sin(K(a;  —  at)).  (75) 

We  solved  this  problem  with  the  same  hrst-order  hyperbolic  advection-diffusions 
equation  given  as  in  Eqs.  (71)  and  (72).  For  each  physical  time,  we  reduced 
the  residuals  by  two  orders  of  magnitude  before  advancing  in  time.  During 
each  time  step,  we  also  relaxed  the  linear  system  using  GS  relaxations  until 
two  orders  of  magnitude  reduction  in  the  linear  system  residuals  was  achieved. 
Shown  in  Fig.  7  are  the  solution  and  solution  gradient  on  =  25  uniform  grid 


Figure  7:  Time-dependent  linear  advection-diffusion  problem(a  =  1,  z/  =  0.03) 
with  periodic  boundary  condition  on  =  25  uniform  nodes.  Time  step  = 
0.01s.  From  left  to  right,  solutions  at  t  =  0.1s,  0.5s,  1.0s. 

for  a  =  1,  1/  =  0.03,  and  k  =  27r.  The  results  (Hyp-ADE)  are  over  plotted 
with  the  exact  solution,  indicating  excellent  accuracy  on  such  a  coarse  grid. 
Similar  results  were  also  obtained  for  larger  grid  systems. 

We  examined  the  convergence  rate  of  this  problem  on  several  grid  systems. 
Given  in  Table  3,  are  the  average  numbers  of  GS  relaxations  per  Newton  itera¬ 
tion  obtained  over  100  time  steps.  Glearly  the  convergence  rate  is  of  O(A^),  not 
0{N‘^)  as  typical  for  numerical  methods  for  the  advection-diffusion  equation. 
Observe  that  for  most  grid  systems  only  one  or  two  Newton  iterations  were 
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Table  3:  Hyperbolic  time-dependent  linear  advection-diffusion  problem  with 
periodic  BC  (a  =  1,  z/  =  0.03)  Average  data  over  100  time  steps  are 
given. (Convergence  criteria:  Residuals  <  10“^) 


Number  of  nodes 

GS  relaxations/Newton  iteration 

Newton  iterations 

25 

4 

3 

100 

11 

3 

300 

33 

2 

500 

55 

2 

1000 

116 

2 

sufficient  to  obtain  accurate  solutions.  We  also  note  that  we  used  At  =  0.01- 
sec  for  all  the  grid  systems;  the  BDF  schemes  are  unconditionally  stable.  The 
time  step  is  orders-of-magnitude  larger  than  that  required  for  conventional 
explicit  schemes,  which  is  limited  by  0{h^).  Of  course,  conventional  implicit 
schemes  also  allow  unconditionally  large  time  steps,  but  it  requires  0{N‘^) 
convergence  in  an  iterative  linear  solver  and  potentially  a  much  larger  num¬ 
ber  of  outer  iterations  as  well  if  the  exact  linearization  is  not  possible  and 
Newton’s  method  cannot  be  constructed.  Hence,  the  method  developed  here 
has  two  major  advantages  over  conventional  methods:  the  exact  linearization 
(Newton’s  method)  and  0{N)  iterative  convergence  in  the  linear  solver.  The 
latter  advantage  can  be  potentially  huge  with  increase  of  the  grid  system  as 
the  speed-up  factor  is  0{N)  and  thus  grows  for  finer  grids.  For  example,  for 
N  =  1000,  our  scheme  is  at  least  four  orders  of  magnitude  faster  than  the 
conventional  schemes,  which  is  remarkable. 

6.3  Unsteady  Nonlinear  Advection-Diffusion 

Consider  the  time-dependent  nonlinear  viscous  Burgers  equation  with  a  time- 
dependent  source  term: 

dtu  + dxf  =  d^iuua;)  +  g{x,t),  a;e(0,l),  (76) 

where  /  =  u^/2,  u  =  u,  and 

g{x,t)  =ut+^iiu^f)cc- {ulf -u^ul^.  (77) 

The  source  term  has  been  generated  by  the  following  function: 

smh.{^y  ioj  /  u) 


j  +  C,  C  >  1,  (78) 
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so  that  it  is  the  exact  solution  to  Eq.  (76).  The  boundary  conditions  are 


u{0,t)  =  C,  (79) 

=  C  +  Ucos{ut),  (80) 

where  uj  is  the  frequency  of  the  oscillation  on  the  right  boundary,  and  U  is  the 
amplitude  of  the  oscillation.  The  constant  C  must  be  greater  than  1  to  keep 
the  diffusion  coefficient  positive. 

We  solved  this  time-dependent  nonlinear  advection-diffusion  equation  with 
the  following  equivalent  hrst-order  hyperbolic  system: 

3  _  ^n-l 

drU  +  4  (u^)  =  d^p  -  —u  + - — - +  g{x,  t),  (81) 

'^drP  =  (d^u-p/u).  (82) 

Shown  in  Fig.  8  are  the  results  obtained  on  iV  =  100  nonuniform  nodes  for 
U  =  1,  CO  =  Ttc/2,  and  C  =  2.  The  scheme  produces  very  accurate  time 


Figure  8:  Time-dependent  nonlinear  advection-diffusion  {a  =  v  =  u)  with 
oscillatory  boundary  condition  u{0,t)  =  2,  u{l,t)  =  2  +  cos(^t)  on  =  100 
nonuniform  nodes.  Time  step  =  0.01s.  From  left  to  right,  solutions  at  t  = 
0.1s,  2.5s,  10.0s. 

evolutions  of  the  solution  [u)  and  the  gradient  {p/io)  on  nonuniform  grids  (see 
also  Fig.  9.) 

The  0{N)  convergence  rate  was  once  again  achieved  for  the  time-dependent 
nonlinear  hyperbolic  advection-diffusion  system.  This  is  given  in  Table  4, 
where  the  average  number  of  iterations  were  obtained  over  1000  time  steps 
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Figure  9:  Time-dependent  nonlinear  advection-diffusion  (a  =  u  =  u)  with 
oscillatory  boundary  condition  u{0,t)  =  2,  u{l,t)  =  2  +  cos(^t)  on  =  100 
nonuniform  nodes.  Time  step  =  0.01s 


(over  17  periods).  Note  also  that  the  method  converged  at  four  Newton  iter¬ 
ations  for  all  grids  considered.  Again,  a  constant  time  step  (At  =  0.01-sec) 
was  used  for  all  the  grid  systems.  And  the  method  is  substantially  more  ef- 
hcient  here  also  than  conventional  explicit /implicit  schemes  as  discussed  for 
the  linear  case  in  Section  6.2. 

Table  4:  Hyperbolic  time-dependent  nonlinear  advection-diffusion  problem 
with  BDF2.  Average  data  over  100  time  steps  are  given.  (Convergence  criteria: 
Residuals  <  10“^) 


Number  of  nodes 

GS  relaxations/Newton  iteration 

Newton  iterations 

100 

88 

4 

300 

263 

4 

500 

440 

4 

1000 

882 

4 

6.4  Accuracy  Verification 

To  demonstrate  the  overall  order  of  accuracy  of  our  time-dependent  hyperbolic 
advection-diffusion  scheme,  we  computed  the  Li  = 

as  Loo  =  —  Ui)  norms  of  errors  on  a  series  of  nonuniform  grids  of 

N  nodes,  and  time  step  increments  At.  For  each  case,  we  terminate  the  solver 
when  the  residuals  were  dropped  by  two  orders  of  magnitude.  We  remark 
that  in  some  practical  applications  two  orders  of  magnitude  reduction  in  the 
residuals  may  not  be  sufficient  to  achieve  the  appropriate  order  of  accuracy. 
But  for  the  cases  considered  here,  terminating  at  even  nine  orders  of  magnitude 
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resulted  in  similar  order  of  accuracy.  The  spatial  accuracy  was  obtained  by 
comparing  several  nonuniform  grid  systems  on  a  fixed  At,  while  the  temporal 
accuracy  was  obtained  by  comparing  various  At  on  a  fixed  nonuniform  grid 
system.  We  also  included  boundary  nodes  in  both  Li  and  L^o  norms.  For 
unsteady  cases,  the  Li  and  L^o  norms  were  evaluated  at  t  =  1.0-sec. 


2  -3 


Figure  10:  Loo  error  obtained  for  time-dependent  linear  hyperbolic  advection- 
diffusion  system  scheme.  Left:  Spatial  accuracy,  Right:  Temporal  accuracy. 

Figures  10  and  11  show  the  L^o  error  convergence  results  for  both  linear 
and  nonlinear  time-dependent  hyperbolic  advection  diffusion  system,  where  h 
is  the  representative  mesh  spacing  defined  by  h  =  /(iV  —  1).  For  discussion 
purposes,  we  only  present  the  accuracy  plots  for  Re  =  1  results;  similar  results 
were  obtained  for  other  Reynolds  numbers.  These  results  show  that  our  scheme 
is  second-order  accurate  for  all  the  variables  and  the  gradients  (note  that  again 
that  our  flux  integral  is  exact)  at  all  the  grid  nodes  including  the  boundary 
nodes  (see  also  Tables  5,  6,  7,  and  8).  We  note  that  it  is  natural  for  the  error 
convergence  in  the  temporal  space  to  deteriorate  when  the  0{AR)  <<  0{h‘^) 
as  well  as  when  At  is  too  large  for  the  time  discretization  to  be  in  asymptotic 
range;  i.e.  a  lower  order  of  accuracy  for  coarse  discretization  in  time  is  typical. 
These  limits  are  also  included  in  the  temporal  accuracy  plot  for  completeness. 

Table  5:  Spatial  accuracy  for  the  linear  advection-diffusion  problem  with  oscil¬ 
latory  boundary  condition  on  several  nonuniform  grid  systems  {At  =  0.001s). 


Number  of  nodes 

Li  error  of  u 

Loo  error  of  u 

Order 

Li  error  of  p 

Loo  error  of  p 

Order 

25 

9.65E-04 

3.44E-03 

6.58E-03 

1.71E-02 

50 

2.28E-05 

8.98E-04 

1.94 

1.08E-03 

4.22E-03 

2.02 

100 

6.10E-05 

2.13E-04 

2.08 

3.08E-04 

9.95E-04 

2.08 

200 

1.59E-06 

5.04E-05 

2.08 

8.76E-05 

2.35E-04 

2.08 

300 

7.58E-06 

2.06E-05 

2.21 

4.87E-05 

9.53E-05 

2.23 

26 


Figure  11:  L^o  error  obtained  for  time-dependent  nonlinear  hyperbolic 
advection-diffusion  system  scheme.  Left:  Spatial  accuracy,  Right:  Temporal 
accuracy. 


Table  6:  Temporal  accuracy  for  the  linear  advection-diffusion  problem  with 
oscillatory  boundary  condition  on  a  nonuniform  grid  (N=100)  with  hMax.  = 


0.045. 


At 

Li  error  of  u 

Loo  error  of  u 

Order 

Li  error  of  p 

Loo  error  of  p 

Order 

0.500 

3.01E-01 

6.63E-01 

2.46E+00 

4.09E+00 

0.100 

2.32E-02 

7.55E-02 

1.35 

1.25E-01 

2.32E-01 

1.78 

0.050 

4.37E-03 

1.72E-02 

2.13 

3.57E-02 

6.28E-02 

1.88 

0.025 

9.66E-04 

3.16E-03 

2.44 

1.16E-02 

2.26E-02 

1.47 

0.010 

1.48E-04 

2.99E-04 

2.57 

2.21E-03 

4.63E-03 

1.73 

0.001 

6.02E-05 

2.13E-04 

0.15 

3.05E-04 

9.95E-04 

0.67 

Table  7:  Spatial  accuracy  for  the  nonlinear  advection-diffusion  problem  with 
oscillatory  boundary  condition  on  several  nonuniform  grid  systems  {At  = 

0.001s). 

Number  of  nodes 

Li  error  of  u 

Loo  error  of  u 

Order 

Li  error  of  p 

Loo  error  of  p 

Order 

10 

25 

1.83E-03 

4.36E-04 

8.41E-03 

1.50E-03 

1.88 

1.20E-02 

2.79E-03 

4.46E-02 

7.29E-03 

1.98 

50 

9.00E-05 

3.74E-04 

2.00 

5.48E-04 

1.78E-03 

2.03 

100 

2.80E-05 

9.39E-05 

1.99 

1.41E-04 

4.36E-04 

2.03 

200 

8.19E-06 

2.33E-05 

2.01 

4.81E-05 

1.07E-04 

2.03 
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Table  8:  Temporal  accuracy  for  the  nonlinear  advection-diffusion  problem 
with  oscillatory  boundary  condition  on  a  nonuniform  grid  system  (N=100) 
with  hMax.  =  0.045. 


At 

Li  error  of  u 

Loo  error  of  u 

Order 

Li  error  of  p 

Loo  error  of  p 

Order 

0.200 

6.46E-02 

1.55E-01 

4.97E-01 

6.83E-01 

0.100 

1.52E-02 

4.67E-02 

1.73 

9.00E-02 

1.54E-01 

2.15 

0.050 

2.03E-03 

7.40E-03 

2.66 

1.99E-02 

3.45E-02 

2.16 

0.020 

3.52E-04 

7.32E-04 

2.52 

5.01E-03 

1.03E-02 

1.32 

0.010 

1.34E-04 

3.20E-04 

1.19 

1.70E-03 

3.06E-03 

1.75 

0.005 

6.68E-05 

1.41E-04 

1.18 

6.20E-04 

9.90E-04 

1.63 

0.001 

2.80E-05 

9.39E-05 

0.25 

1.41E-04 

4.36E-04 

0.51 

7  Conclusions 

In  this  paper,  we  have  extended  the  hrst-order  hyperbolic  system  method  to 
time-accurate  computations.  The  time  integration  scheme  is  implicit  with  the 
second-order  backward-difference  formula,  and  the  resulting  system  of  time- 
dependent  residual  equations  is  solved  by  a  steady  solver  developed  based  on 
the  hyperbolic  method.  The  steady  solver  is  Newton’s  method  constructed 
from  a  compact  upwind  residual-distribution  scheme  for  linear  and  nonlinear 
advection-diffusion  equations.  The  construction  of  the  numerical  scheme  was 
greatly  simplihed  for  the  nonlinear  equation  by  the  preconditioned  conservative 
formulation.  We  have  demonstrated  that  the  developed  time-accurate  schemes 
are  capable,  at  every  physical  time  step,  of  producing  second-order  accurate 
solution  and  gradient  on  highly  rehned  nonuniform  grids  by  less  than  hve 
Newton  iterations  with  the  linear  convergence  of  the  block  GS  relaxations  for 
all  test  cases  considered.  Thus,  the  notable  features  of  the  hyperbolic  method 
have  been  shown  to  carry  over  to  time-dependent  problems. 

The  developed  methodology  based  on  the  second-order  backward  Euler 
time  integration  can  be  extended  to  practical  one-dimensional  problems  [7-11], 
bringing  significant  improvements  in  both  accuracy  and  efficiency.  Extension 
to  higher-order  accuracy  is  also  possible  by  improving  the  residual  evaluation 
and  higher-order  backward-difference  formulas.  It  is  relatively  simple  in  one 
dimension  because  the  integration  of  the  flux  balance  term  is  exact  in  one 
dimension  as  mentioned  in  Section  3.1.  In  higher  dimensions,  although  the 
construction  of  time-accurate  schemes  is  as  simple  as  presented  in  the  paper, 
it  is  known  that  a  high-order  quadrature  would  be  required  for  the  flux  balance 
integration  in  order  to  achieve  the  equal  order  of  spatial  accuracy  in  the  solu¬ 
tion  and  the  gradient  [2].  Although  not  considered  in  the  present  paper,  the 
construction  of  non-oscillatory  schemes  is  an  important  area  of  development 
for  problems  with  discontinuities  such  as  shock  waves. 
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